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Abstract: We analyze an M/M/l queue with a service discipline in which customers, upon arriving when 
the server is busy, search a sequence of stations for a vacant station at which to wait, and in which the 
server, upon becoming free when one or more customers are waiting, searches the stations in the same order 
for a station occupied by a customer to serve. We show how to find complete asymptotic expansions for all 
the moments of the waiting time in the heavy traffic limit. We show in particular that the variance of the 
waiting time for this discipline is more similar to that of last-come-first-served (which has a pole of order 
three as the arrival rate approaches the service rate) than that of first-come-first-served (which has pole of 
order two). 
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1. Introduction 

We consider the M/M/l queue (with independent exponentially distributed interarival times, indepen- 
dent exponentially distributed service times, and a single server) with various service disciplines. We shall 
be interested mainly in the "heavy traffic" limit, A — » 1, where A is the arrival rate (measured in units of the 
service rate), and all asymptotic statements in this paper refer to this limit. 

It is well known (see Little [L]) that the average waiting time Ex[VF] (W is the length of interval from 
arrival to commencement of service) does not depend on the service discipline (the rule used to determine 
which waiting customer is served next when the server becomes free). We have 

The variance of W, however, (and more generally its higher moments) does depend on the service 
discipline. Kingman [K] has shown that, among all service disciplines, "first-come-first-served" (FCFS) 
minimizes the variance of W . We have 

Var[T4Ws] = ~ (1.2) 

(see for example Riordan [R2, pp. 102-103]). Tambouratzis [T] has shown that "last-come-first-served" 
(LCFS) maximizes this variance. We have 



Vl »-[""^ = H2 (i-x/ ] ~JT^W (L;J,) 



(see for example Riordan [R2, pp. 106-109]; LCFS was first analyzed by Vaulot [V2]). We note that the 
difference between FCFS and LCFS is qualitative, in that Var[WLCFs] ^ as a P°l e 01 or der three at A = 1, 
whereas Var[WpCFs] has a pole of order only two there. 

Another service discipline that has been studied is "random-order-of-service" (ROS), first successfully 
analyzed by Vaulot [VI]. One of the motivations for studying ROS was stated by Riordan [Rl]: 

"In many switching systems it is not feasible to fully realize this ethical ideal of first come, 
first served, and it has long been of interest to determine delays on another basis. The contrasting 
assumption is of calls picked at random, which is again an idealization but in large offices appears 
to be called for, as a bound for the service actually given." 

This statement suggests that (1) in practical systems of that era (around 1953) it was not possible to 
keep track of the order of arrival, (2) ROS was analyzed as a substitute for the service discipline actually 
implemented, and (3) it was hoped that the performance of ROS would approximate that of the discipline 
actually implemented. For ROS, we have 

(see for example Riordan [R2, pp. 103-106]). Comparing (1.4) with (1.2) and (1.3), we see that ROS is 
qualitatively "more like" FCFS than LCFS, in that Var[VFRos] has a pole of order two rather than three, 
though its coefficient is larger than that of FCFS by a factor of three. 
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In the early 1950s, telephone switching systems were electromechanical and did not employ random- 
ization beyond that present in the arrival and service processes. In this paper we shall analyze a service 
discipline that we shall call "fixed-order-of-search" (FOS). For this discipline, there is an infinite sequence 
Wi, W2, ... of "waiting stations", each of which can be either "vacant" or "occupied". A customer arriving 
when the server is busy searches these stations in increasing order of their indices and occupies the first 
vacant station it finds. When the server becomes free and one or more customers are waiting, it searches 
the stations in the same order and serves the customer waiting at the hrst occupied station it finds, thereby 
vacating that station. This discipline was introduced by Eschenfeldt, Gross and Pippenger [E]. Like FCFS 
and LCFS (and unlike ROS), it does not employ randomization beyond that present in the arrival and service 
processes. If the server were to search the stations in the reverse order to that used by arriving customers 
(serving the customer waiting at the occupied station with the largest, rather than the smallest, index), the 
result would be LCFS, with its concomitant maximum variance for the waiting time. The choice of the same 
order of search for both customers and the server thus represents an attempt to improve upon LCFS, while 
still using a fixed order of search in each case. 

We shall give an exact formula for Var[WFos] : 

where ^q\x) is the g-trigamma function (defined below). We shall indicate how similar, but increasingly 
more complicated, formulas can be derived for the higher moments of Wfos- We shall also give an asymptotic 
formula for Var[WFOs]: 

Var[I4W] ~ (1-6) 

where ((s) = X)n>i l/ nS i s the Riemann zeta function and £(2) = 7r 2 /6 (see for example Whittakcr and 
Watson [W2, pp. 265-280]). We shall also indicate how complete asymptotic expansions (with error terms of 
the form 0((1 — X) R ) for any R) can be derived for the variance of Wfos, as well as for the higher moments. 
Comparing (1.6) with (1.2) and (1.3), we see that FOS is qualitatively "more like" LCFS than FCFS, in that 
Var[WFos] has a pole of order three rather than two, though its coefficient is smaller than that of LCFS by 
a factor of 4 - 2 C(2) = 4 - tt 2 /3 = 0.7101 .... 

Eschenfeldt, Gross and Pippenger [E] initiated the study of FOS, determining the distribution of the 
index I of the station Wj at which a newly arriving customer waits (where I = if the server is idle at the 
time of the arrival): we have Pr[J > 0] = 1 and 

Pr[/ > .1 = (1.7) 
for i > 1. The moments of / can be expressed in terms of the sums 

i>l 

for which we have the asymptotic formulas 

H)(A)~ r ^log T ^ (1.9) 
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and 



Ti(X) 



(1 - A) i+1 



(1.10) 



for I > 1, and the exact formulas 



T (A) = 



V> A (1) + log(l - A) 
log A 



(1.11) 



and 



Ti(X) 




log l+1 X 



(1.12) 



for i > 1. Here ip q (x) = dlogT q (x)/dx is the q-digamma function, the logarithmic derivative of the g-gamma 
function T q (x) = (1 — q) 1 ^ x Y\ n>0 ((l — <7™ +1 )/(l — <7" +a: )) (sec for example Gasper and Rahman [G, p. 16]), 
and ip q l \x) = d l ip q (x) / dx l is the Z-th g-polygamma function. (The sums T";(A), which are called Lambert 
series (see for example Hardy and Wright [H, p. 257]), are the generating functions TJ(A) = X)«>i a i( n ) A n 
for the sums er;(n) = J2d\n ^ l °^ tnc ^" tn P owers of the divisors of of n (see for example Hardy and Wright 
[H, p. 239]).) In terms of the T ; (A), we have 



In Section 2 we shall determine the moment generating M\y (s) function for Wfos (which in what follows 
we shall denote simply W). In Section 3, we shall derive the exact formula (1.5) and the asymptotic formula 
(1.6). We shall also indicate how similar exact and asymptotic formulas can be found for the higher moments 
of W . Finally, we shall indicate how these asymptotic formulas can be extended to complete asymptotic 



2. The Generating Functions 

Consider the random process whose state variable J denotes the number of customers in the system. 
The random variable J is zero during an idle period (interval of time during which the server is idle). It 
is incremented whenever a customer arrives, and decremented whenever a customer departs (that is, at the 
termination of a service interval). Arrivals occur in a Poisson process with rate A. During a busy period 
(interval of time during which J > 1), departures occur in an independent Poisson process at rate 1. Thus, 
during a busy period, "transitions" , by which we mean arrivals and departures together, occur in a Poisson 
process at rate 1 + A. Furthermore, during a busy period, the probability that the next transition will be 
an arrival is p = A/(l + A), and the probability that it will be a departure is q = 1/(1 + A). In this section 
we shall study the distribution of the random variable N, defined as the number of transitions that occur 
between the arrival of a customer (excluded) and the departure that initiates its service interval (included). 
Specifically, we shall determine the probability generating function g(t) = ^2 n>0 Pr{N — n] t n for N. We 
have N = if the arrival initiates a busy period, and N > 1 if it occurs during a busy period. 

The index I of the station Wi at which a newly arriving customer waits has the distribution given by 
(1.7). As a first step to determining g(t), we shall determine the conditional generating function gi(t) = 
J2 n >a — n\ I — i]t n . We begin with two special cases. If i = 0, the arrival initiates a busy period, so 
N = and go(t) = 1. If i = 1, the customer waits at W\ and will be served as soon as the next departure 
occurs. The number of transitions preceding and including this next departure has a geometric distribution, 




(-1) 



m—l—l 




expansions (with error terms of the form 0((1 — X) R ) for any R) for these quantities. 



and gi (t) =E„>i:P"~V" = qt/(l-pt). 
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For the general case, we consider the random process whose state variable K denotes the number of 
vacant stations among Wi, . . . , W,. Since the customer in question waits at station Wj, we have K = 
immediately after the arrival of that customer. Furthermore, the first time thereafter at which K = i 
coincides with the beginning of the service interval for that customer, and thus occurs after exactly N 
transitions have occurred. When 1 < K < i — 1, the random variable K is incremented by a departure 
(because the next customer served will be waiting at one of the stations under consideration) and decremented 
by an arrival (because the arriving customer will wait at one of these stations). If however K — 0, a departure 
will increment K, but an arrival will leave K unchanged (because it will have to wait at a station beyond 
Wi). Thus the process determining N given / = i is a random walk with one "reflecting barrier" (at K = 0) 
and one "absorbing barrier" (at K = i), as shown in Figure 1. 

p 



siari p 




Figure 1. The process determining N, given that I = i. 
Non-terminal states are labeled with values of the random variable K. 

The number N of steps to absorption in this process has been studied by Weesakul [Wl], who shows 
that the generating function is given by 

„. (t) = iim-m) 

yt{ ' Q(ty+i - p(t)<+i - P t(Q(ty - p(ty) ' K ' ' 

where 



Pit) = WlUgg 



and 



Q(t) = i- + 



2 

We note that for t = 1 we have P(l) = p, Q(l) = q and gi(l) = 1. 

We can now express the unconditional generating function g(t) by using summation by parts: 



ff(t) = 5>(*) Pr[/ = i] 

i>0 

= Y,9i(t) (Pr[J>i]-Pr[J>* + l]) 

= go(t) Pr[J > 0] + - ft-i(t)) Pr[J > i] 

i>l 

= 1 + S(ft(*)-ft-l(*)) Pr[/> 4 ], 



t>l 



because go(t) = Pr[7 > 0] = 1. Thus, using (1.7), we have 

(Vff)i(t) A 



ff(t ) = ! + (!- A) S™^, (2-2) 
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where (Vg)i(t) = gt(t) — gi-i(t) denotes the backward difference of g%{t). We note that for t = 1 we have 
(V ff )i(l)=0, so = 1. 

We are now ready to drive the moment generating function M\y(s) — Ex[e sW ] for W. Each intertran- 
sition time X is exponentially distributed with mean 1/(1 + A), so the moment generating function for X is 
Mx(s) = (1 + A)/(l + A — s). The moment generating function for the sum ^2 1<k<n Xk of n independent 
intertransition times X\, . . . ,X n is Mx(s) n = ((1 + A)/(l + A — s)) n . Thus the waiting time W, which is 
the sum of the random number N of independent intertransition times has the moment generating function 

M w (s) = Pr[^ = n] M x {s) n 

n>0 

= g(M x (s)) (2.3) 

= 9 {iTxh)- 



3. The Moments 

In this section we shall derive the mean and variance for N and W, and indicate how to derive the 
higher moments as well. 

For the mean of N, we use the formula Ex[iV] = <?'(1). We have 



(l + A) (i(l-A)-A(l-A')) 
(1 - A)^ 



so 



Thus we have 



(v^(D = ™^- 



EX[7V] = g'(l) 



(1_A) 4- l-A* 



»>i 



l-A' 



_,, (l + A) (l-A*) A 

= (l + A) A 
l-A ' 



(3.1) 



For the mean of W, we use the formula for the expectation of the sum of a random number N of indepen- 



dent, identically distributed random variables X,X\,X2, ■ ■ ■ ■ Ex 
intertransition time X satisfies Ex[X] = 1/(1 + A), we have 

Ex[W] = Ex[AT] Ex[X] 
(l + A) A 1 



Ei<k<N X k = Ex[iV]Ex[X]. Since the 



l-A l+A 
A 



l-A' 



in accordance with (1.1). 
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For the variance of N, we begin by using the formula for the second factorial moment: Ex[N(N — 1)] 
g"(l). We have 

„ _ (1 - A) 2 (l + A) 2 » 2 - (1 - A)(l + A)(l - 10A - 3A 2 ) % 
9l ( ' ~ (1 - A) 4 

(6A(1 - A)(l + A) 2 i + 2A(1 + A)(l + 4A - A 2 ) - 2A 2 (1 + A) 2 A*)(l - A*) 

so 

(2(1 - A)(l + A) 3 A* + 6»(1 - A) 2 (l + A) 2 - 2(1 - A) 3 (l + A)) (1 - A*) - 4*(1 - A) 2 (l + A) 2 



(Vff")i(l) 



(1-A) 4 
Thus we have 

v ; i>\ »>1 i>l i>l 

2A 2 (1 + A) 2 6A(1 + A) 2 2A(1 + A) 4(1 + A) 2 

(1-A) 3 (1-A) 3 1-A 1-A U >' 



where we have used the definition (1.8) to evaluate the last sum. It follows that 

Var[7V] = Ex[iV(iV - 1)] + Ex[N] - Ex[7V] 2 

_ 2A 2 (1 + A) 2 6A(1 + A) 2 A(l + A) A 2 (l + A) 2 4(1 + A) 2 
(1-A) 3 + (1-A) 3 ~~ 1-A ~~ (1-A) 2 ~~ 1-A 
A(l + A) 2 (6 + A + A 2 ) A(l + A) 4(1 + A) 2 



(1-A) 3 1-A 1-A 



71(A), 



where we have used (3.1), then combined the first, second and fourth terms. For the variance of W, we 
use the formula for the variance of a random number N of independent, identically distributed random 
variables X,X\, X2, Var J2 1<k<N Xk = Var[7V] Ex[Jf] 2 + Ex[7V] Var[X]. Since the intertransition time 
X satisfies Ex[X] = 1/(1 + A) and Var[X] = 1/(1 + A) 2 , we have 

Var[W] = Var[iV] Ex[X} 2 + Ex[N] Vai[X] 

'A(l + A) 2 (6 + A + A 2 ) A(l + A) 4(1 + A) 2 \ 1 A(l + A) 1 



(1-A) 3 1-A 1-A ' 1K "'J (1 + A) 2 ' 1-A (1 + A) 5 

A(6 + A + A 2 ) 4 



(1-A) 3 1-A 
Evaluating Ti(A) using (1.12) yields 



11(A). 



Var W ^( 6 n +A + A2 )-/^W , 
(1-A) 3 (1-A) log 2 A' 

confirming (1.5), whereas using (1.10) yields 

v» rwi A (6 + A + A 2 ) 4 C(2) 

Var[T^] ~ (1 _ A)3 - JT - W , 
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confirming (1.6). 

It is straightforward to generalize the derivations of the mean and variance of W given above to the 
higher moments. We begin by indicating how to derive the higher factorial moments of N. We have 

Ex[iV(7V - 1) • • • (N - m + 1)] = 5 (m) (1) 

" h ' 

1>1 

After differentiating g^t) with respect to t (m times), then evaluating the result at t — 1, and finally 
differencing with respect to i, the result is a bivariate polynomial V(i, u) (with coefficients that are rational 
functions of A) in the variables i and u — A 4 . Dividing this polynomial by 1 — u = 1 — A*, we obtain 
V(i, u) — Q(i, u)(l — u) + TZ(i), with quotient Q(i, u) and remainder lZ(i). We then have 



TZ(i) A 4 

i>i »>i 



Ex[7V(iV - 1) • • ■ (N - m + 1)] = (1 - A) £ Q(i, A 4 ) A 4 + (1 - A) £ y^-. (3.2) 



The first sum in (3.2) can be expressed as a linear combination (with coefficients that are rational functions 
of A) of sums of the form 



»>i 

These sums arc themselves rational functions of A: 



(l-A*)^ 1 ' 



where Ai{x) — J2o<k<i a (l > ^) x h is the Z-th Eulerian polynomial and the a(l,k) are the Eulerian numbers, 
with generating function J2i k>o ^) zk V 1 1^- = z (l — z)/(e v ^~ z ^ - z) (see Comtet [C, p. 245]). The second 
sum in (3.2) can be expressed as a linear combination (again with coefficients that are rational functions 
of A) of the sums TJ(A) given by (1.8), with asymptotic formulas given by (1.9) and (1.10), and with exact 
formulas given by (1.11) and (1.12). 

We are now ready to obtain the moments of W. Differentiating the identity (2.3) m times, we obtain 
M&Hs) = F m (M£\s), . . .,M^\s);g^(M x (s)), . . . , g^(M x (s))) , 
where F m (x\, . . . , x m ; yi, . . . , y m ) is the polynomial 

F m (x 1 ,...,x m ;y 1 ,...,y m ) = m! x l JJ (|y) , 

l<l<m ei ,C2 , . . . ,e m l<k<m 

and the inner sum is over all e\, ei-, . . . , e m such that e\ + 62 + ■ ■ ■ + e m = I and e\ + 2 C2 + ■ ■ ■ + m e rn = m 
(see for example Comtet [C, p. 137]). Evaluating at s = and using M x (0) = 1 and M x \o) = 1/(1 + A) ( 
for I > 1, we have 

M&\0) = F m (gW(l),...,gW(l); T ^,. 



1 



(1 + A)* 



= nTxr S E ( eie2 m e ) 

»- 1 + A J l<i<m ei, e2 ,...,e m V e l > e 2 , • • • , W 

v ; KKm 
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where the { ™ } are the Stirling numbers of the second kind, with the generating function 
Em>;>0 { 7}y l z " 1 / 11 = e v(e '- 1) (see for example Comtet [C, pp. 206-207]). Thus we have 

Ex[VF m ] =M^ l) (0) 

V ; l<l<m 

where the g( l \l) are the factorial moments of N determined in the preceding paragraph. 

The asymptotic formulas given above for the moments of W can be extended to complete asymptotic 
expansions, with error terms of the form 0((1 — X) R ) for any R. Any rational function of A has a Laurent 
series around A = 1, which will serve an an asymptotic expansion as A — > 1 as well. Thus the only 
remaining problem is to find asymptotic expansions for the sums I] (A). These expansions have been given 
by Eschenfeldt, Gross and Pippenger [E]. We have 

T 0( A)~-iog— — (r+1)(r+1)! — - 

where 7 = 0.5772... is Eulcr's constant, B r is the r-th Bernoulli number, defined by t/(e t — 1) = 
E r >o^ r *V r ' ( see f° r example Roman [R3, p. 94], and h — — log A has the expansion 

fc=-log(l-(l-A)) 

r ' 

r>l 



so that its reciprocal has the expansion 



1 1 v (-iy c r (i- xy 

h l-A 4r( 



r 

r>0 



where C r is the r-th Bernoulli number of the second kind, defined by tj log(l + t) = J2r>o t r /r\ (see for 
example Roman [R3, p. 116]). (These numbers are also called the Cauchy numbers of the first kind, and are 
given by C r = x(x — 1) • • • (x — r + 1) dx; see for example Comtet [C, pp. 293-294].) For I > 1, we have 

^ /x x (-l) r+ ' _1 B r B r+ i h r ~ l 

T ^ - 2-, h(7To • 

We note that, if I is odd, then this expansion has only finitely many terms (because B r = for odd r > 3). 
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